function [FitParameters, StandardErr, RedChiSqr, exitflag] = ...
    lsq_fixedp(Function, InitialParameters, FixedParameters, LowerLimits, UpperLimits)
%[FitParameters, StandardErr, RedChiSqr] = ...
%    lsq_fixedp(Function, InitialParameters, FixedParameters, LowerLimits, UpperLimits)
% For easy use of lsqnonlin with some parameters fixed. Like lsqnonlin,
% this accepts a function, where the square fo this function is going to be
% minimized. So, e.g., if you want to fit points {y_i} to a model function
% f(P, x) for points {x_i}, you would need to supply function g(P) =
% [f(P,x_1) - y_1, f(P, x_2) - y_2, ...]
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Check arguments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~ischar(Function) && ~isa(Function, 'function_handle')
    errorStruct.message = 'Argument Functions was not a string of characters function handle.';
    errorStruct.identifier = 'lsq_fixedp:ArgType';
    error(errorStruct);
end

if ~exist('FixedParameters', 'var') || isempty(FixedParameters)
    FixedParameters = zeros(1,length(InitialParameters));
end

if length(FixedParameters) ~= length(InitialParameters)
    FixedParameters = zeros(1, length(InitialParameters));
    warning('FixedParameters was a different size than InitialParameters. Set to zeros.');
end

if any(FixedParameters ~= 0 & FixedParameters ~= 1)
    FixedParameters(FixedParameters ~= 0 & FixedParameters ~= 1) = 0;
    warning('One or more fixed parameters was not a zero or one. Set offending parameters to zero');
end

if ~exist('UpperLimits', 'var') || isempty(UpperLimits)
    UpperLimits = inf(length(InitialParameters),1);
end

if length(UpperLimits) ~= length(InitialParameters)
    UpperLimits = inf(length(InitialParameters),1);
end

if ~exist('LowerLimits', 'var') || isempty(LowerLimits)
    LowerLimits = -inf(length(InitialParameters),1);
end

if length(LowerLimits) ~= length(InitialParameters)
    LowerLimits = -inf(length(InitialParameters), 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ischar(Function)
    fn_handle = str2func(Function);
elseif isa(Function, 'function_handle')
    fn_handle = Function;
else
    errorStruct.message = 'Function was not a function handle or a string';
    errorStruct.identifier = 'lsq_fixedp:value';
    error(errorStruct);
end

%Now create a new anonymous function that doesn't include any of our fixed
%parameters. This is all we need for lsqnonlin
FullFunctionHandle = @(Pred) fn_handle(Reduced2FullParams(Pred,InitialParameters,FixedParameters));

%Do fitting and extract parameters to return.
ReducedInitialParameters = Full2ReducedParams(InitialParameters,FixedParameters); %define this here because I need to know how many there are.
redlowerbs = Full2ReducedParams(LowerLimits,FixedParameters);
redupperbs = Full2ReducedParams(UpperLimits,FixedParameters);
    %use lsqnonlin function
    [RedFitParameters, resnorm, residual, exitflag, output, lambda, jacobian] = ...
        lsqnonlin(FullFunctionHandle, ReducedInitialParameters, redlowerbs, redupperbs);
    RedChiSqr = resnorm / (length(FullFunctionHandle(RedFitParameters)) - length(RedFitParameters));
    %get standard error from confidence intervals
%     NinetyFivePer_ConfInterval = nlparci(ReducedInitialParameters, residual, 'jacobian', jacobian);
    NinetyFivePer_ConfInterval = nlparci(RedFitParameters, residual, 'jacobian', jacobian);
    StandardErr = (NinetyFivePer_ConfInterval(:,2) - NinetyFivePer_ConfInterval(:,1)) / 3.92;
    %convert back to original params
    FitParameters = Reduced2FullParams(RedFitParameters, InitialParameters, FixedParameters);
    StandardErr = Reduced2FullParams(StandardErr, InitialParameters, FixedParameters);

end

%Helper functions to deal with fixed parameters
function [FullParams] = Reduced2FullParams(RedP,InitParams,FixedParams)
FullParams = zeros(size(InitParams));
FullParams(FixedParams == 1) = InitParams(FixedParams == 1);
FullParams(FixedParams == 0) = RedP;
end

function [RedParams] = Full2ReducedParams(FullParams,FixedParams)
RedParams = FullParams(FixedParams == 0);
end